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Abstract 



We consider the general spatial three body problem and study the 
dynamics of planetary systems consisting of a star and two planets 
which evolve into 2/1 mean motion resonance and into inclined orbits. 
Our study is focused on the periodic orbits of the system given in a 
suitable rotating frame. The stability of periodic orbits characterize 
the evolution of any planetary system with initial conditions in their 
vicinity. Stable periodic orbits are associated with long term regular 
evolution, while unstable periodic orbits are surrounded by regions 
of chaotic motion. We compute many families of symmetric periodic 
orbits by applying two schemes of analytical continuation. In the 
first scheme, we start from the 2/1 (or 1/2) resonant periodic orbits 
of the restricted problem and in the second scheme, we start from 
vertical critical periodic orbits of the general planar problem. Most of 
the periodic orbits are unstable, but many stable periodic orbits have 
been, also, found with mutual inclination up to 50° - 60°, which may 
be related with the existence of real planetary systems. 

keywords 2/1 resonance, 3D general three body problem, periodic orbits, 
vertical stability, planetary systems. 
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1 Introduction 



The dynamics of planetary systems consisting of two massive planets can be 
studied by considering the general three body problem (GTBP). Many of 
such systems seem to be locked in mean motion resonances and particularly 
in 2/1 resonance (e.g. GJ 876, HD 82943, HD 73526 and 47 Uma). Such 
resonances are related with periodic orbits of the three body problem in a 
rotating frame and are very important, since they are associated with regions 



of stability and instability in phase space ( 


Psychoyos and Hadjidemetriou 


2005 


Gozdziewski et al. 


2005 


Voyatzis 


2008 


). Particularly, the phase space 



structure near the 2/1 resonance and its connection with periodic evolution 
is studied in Michtchenko et al. ( 2008a|[b ) and Michtchenko and Ferraz-Mello 
(2011). Also, it has been shown that stable periodic orbits can drive the 



migration process of planets 


Lee and Peale 


2002 


Ferraz-Mello et al. 


2003 


Hadjidemetriou and Voyatzis 


2010 


)• 



In the previous years, most of the studies of planetary evolution in mean 
motion resonances were restricted in planar motion. In such cases, reso- 
nant periodic solutions can be found in various ways, see e.g. 
Haghighi pour et al. 2003; [Voyatzis and Hadjidemetriou 2005 



^Varadi 


1999 


Michtchenko 



et al.^2006a ). It is shown that periodic orbits are classified into various config- 



urations, where some of them are stable and others unstable. Setting initial 
conditions close to a stable periodic orbit the planetary system evolves reg- 
ularly, since in phase space and around stable periodic orbits invariant tori 
exist. Along this motion the resonant angles show librations around the val- 
ues that correspond to the particular periodic orbit (see e.g. Voyatzis (2008) 
and Michtchenko and Ferraz-Mello (2011)). 

Concerning the dynamics of multiple planetary systems in space, mu- 
tual inclinations are a very interesting feature. The dynamics of such sys- 
tems can be studied, for instance, by examining the phase space structure 
( [Michtchenko et al.||2006b| or the Hill's stability ( jVeras and Armitage|[2004l ). 
Numerical simulations show that inclined planetary systems can be stable 



even for relatively high mutual inclinations at systems as 47 Uma (Laughlin 



et al. 2002) and v And (Barnes et al. 2011). Stability of inclined systems 



may be associated with the Kozai resonance (Libert and Tsiganis 2009a) 
or the Lagrangian equilibrium points (Schwarz et al. 2012). Some possible 
mechanisms that lead to the excitation of planetary inclinations are the plan- 



etary scattering (Marzari and Weidenschilling 2002 



the differential migration (Thommes and Lissauer 



nis 



Chatterjee et al. 2008), 



2003 Libert and Tsiga- 



2009b Lee and Thommes 2009) and the tidal evolution (Correia et al. 



2011). In these mechanisms, resonance trapping is possible to occur and the 



evolution of the system may be associated with particular inclined periodic 
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orbits. 

In this paper, we study resonant periodic orbits of the spatial three body 
problem that model a possible planetary system in 2/1 resonance. We at- 
tempt to determine all possible periodic configurations and especially the sta- 
ble ones. We apply the method of analytic continuation described and proved 
in ( Ichtiaroglou and Michalodimitrakis|1980 ). Based on this method, families 
of periodic orbits of the general three body problem in space have been com- 



puted by 


Michalodimitrakis 


(1979b) 


Michalodimitrakis 


(1980 


), 


Michalodim- 


itrakis 


(' 


1981 


) and 


Katopodis et al. 


(1980 


). However, these families are not 



associated with planetary resonant dynamics and stability aspects are not 
sufficiently discussed. 

In Section 2, we present the particular TBP model and refer to periodic 
conditions, stability and resonances. In Section 3, we demonstrate typical 
examples of periodic orbits, we present the evolution of orbits starting from 
the vicinity of these periodic orbits and discuss their long-term stability. In 
section 4, we compute and present families of resonant periodic orbits and 
indicate their stability. Some miscellaneous cases are discussed in section 5 
and we conclude in section 6. 



2 The general TBP in space and periodic or- 
bits 

2.1 The model 

Assume a planetary system consisting of a star 5* and two planets Pi and P2 
with masses mo, mi and 1712, respectively, moving in the space under their 
mutual gravitational attraction. We consider these bodies as point masses 
and let O be their center of mass. The three bodies are isolated, thus their 
center of mass is fixed with respect to any inertial system of reference and its 
angular momentum vector, L , is constant. We now introduce a particular 
inertial coordinate system, OXYZ, such that its origin coincides with their 
center of mass O and its Z-axis is parallel to L. The system is described by 
six degrees of freedom defined for instance, by the position vectors of the two 
planets. 

We can reduce the number of degrees of freedom by introducing a suit- 



able rotating frame of reference. Following Michalodimitrakis (1979a), we 
introduce the rotating frame Gxyz, such that: 

1. Its origin coincides with the center of mass G of the bodies S and Pi. 

2. Its 2;-axis is always parallel to the Z-axis. 
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3. 5* and Pi move always on a;z-plane. 

We define Xi, yi, Zi {i = 1, 2, 3) as the Cartesian coordinates of tlie bodies 
witli respect to Gxyz and an angle 9 as the angle between X and x axes 
assuming always that 6{0) = 0. According to the definition of Gxyz, it 
always holds that yi = 0. Considering the normalized gravitational constant 
G = 1 and setting m = mo + mi + m2 (thus 27i time units (t.u.) correspond 
to one year), the Lagrangian of the system in the rotating frame of reference 



is given by Michalodimitrakis (1979a): 



£ = ^(mo + mi)[a{xl + z\ + x\e'^) + h\{x\ + yl + zj) 
+9\xl + yl) + 29{x2m - X2y2)]] - V, 
where V = a = mi/mo, b = m2/m, 

rli = {l + a)\xl + zl), 

= {axi + X2f + yl + {azi + Z2Y , 



(1) 



12 



(Xi - X2f +yl + {Zi - Z2f. 



The angle 6 is ignorable and, consequently, the angular momentum pg = 
82/ 86 is constant and given by 

Pe = ("mo + ■mi)[axl6 + b[9{xl + yl) + {x2y2 - ^2^/2)]] = const. (2) 

Solving the Eq. ^ with respect to 6 we find 



Pe 



- h{x2y2 - X2y2) 



ax^ + 6(X2 + t/a) 

As far as the components of the angular momentum vector are concerned, 
due to the choice of the inertial system it always holds: 



Lx = (fTio + mi)[h{y2Z2 - 2/22:2) - ^(axi^i + hx2Z2)] = 0, 

Ly = (mo + mi)[a(i;i2;i - xiii) + b{x2Z2 - a;2i2) - b9y2Z2] = 



(4) 



and 



Lz = {mo + mi)[6(x2?/2 - X2y2) + 0[axl + b{xl + yl)]], 
which coincides with integral ([2]). The constrains defined by Eqs. ^ give: 



_ mom2 I y2Z2-y2Z2 X2Z2 

1 rrtim I x-^O xi 



(5) 

Zl = ^^,ix2Z2-X2Z2-9y2Z2) + f^Z,. 
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Thus, for a particular value of the angular momentum, p^, the position of 
the system is defined by the coordinate Xi of the planet Pi and the posi- 
tion {x2, 1/2,22) of the planet P2. Namely, the system has been reduced to 
four degrees of freedom. We can derive the equations of motion from the 
Lagrangian ([T|: 

^ mom2(xi-X2) mom2{axi+X2) 

(mo+mi)[(a;i-X2)2+j/2+(2i-2^2)^]^''^ {mo+mi)[{axi+X2y+y2 + {azi+Z2)^]^^^ , . 

(0) 

(mo+mi)xi _i_ ^ a2 
[(l+a)2(x2+22)]3/2 J-lP 



mmi(xi— X2) mm,o{axi+X2) 

(mo+mi)[(a::i--a;2)2+J/| + (2i-Z2)2]3/2 („io+m,i)[{aa;i+X2)2+i/|+(azi+22)2]^/2 

+X2e^ + 2y2e + y2e 

(n 

+7/2^2 - 2X2^ - X2e 



^2 (m,)+mj)[(a;-^_a:^2)2+j/|+{2i-Z2)2]^/2 (mo+mi)[(axi+a;2)2+j/| + (a2;i+2;2)2]^/2 



(7) 



(8) 



" mmi{zi-Z2) mmo{axi+X2) /q\ 

2 ~ {jTto+mi)[(a;i-X2)2+?/^+(2i-2;2)2]3/2 (^mo+mi)[{axi+X2)^+y^ + {azi+Z2)^]'^/'^ ^'^^ 

We should remark some essential differences between the restricted and the 
general 3D-TBP described in the rotating frame. In the general problem, 
the primaries do not remain fixed on the rotating axis Gx, but move on the 
rotating plane xz. Also, the rotating frame Gxyz does not revolve with a 
constant angular velocity, see Eq. ([3]), and its origin G does not remain fixed 
with respect to the inertial system OXY Z. 

In our numerical computations we should always take that Pi is the inner 
planet, P2 is the outer one and the total mass of the system is m = 1. Also, 
we should always take that one of the planets (either the inner or the outer) 
has the mass of Jupiter, namely mi = 10~^, with either i = 1 or i = 2. 



2.2 Periodic orbits 

Considering the rotating frame and the Poincare map y2 = 0, periodic orbits 
are defined as fixed or periodic points of the Poincare map satisfying the 
conditions 

a;i(0) = xi(T), X2(0) = X2(T), ^2(0) = Z2{T), 

xi(0) = Xi(T), X2(0) = X2{T), i2(0) = i2(T), 2/2(0) = UT), ^ ' 

provided that 7/2(0) = y2{T) and T is the period. Due to the existence of the 
Jacobi integral, one of the above conditions is always fulfilled, when the rest 
ones are fulfilled. 
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There exist four transformations under which the Lagrangian ([T]) remains 
invariant and, subsequently, there exist four symmetries for the orbits of 
the 3D-GTBP in the rotating frame (Michalodimitrakis 1979a[ ). A periodic 



orbit defined by (10) is called symmetric, if it is invariant under one of the 
symmetries of the system. In the planetary TBP, there exist periodic orbits 
which are symmetric with respect either to the x2;-plane or to the x-axis. 

A periodic orbit obeying the symmetry with respect to x2;-plane has two 
perpendicular crossings with the xz-plane. Planet, Pi, moves on that plane 
and, without loss of generality, we may consider that planet, P2, is also on 
the a;2;-plane at t = 0. Thus, the initial conditions of such a periodic orbit 
should be 

xi(0) = a;io, X2(0) = X20, 2/2(0) = ^2(0) = ^20, /..n 

i;i(0) = 0, X2(0)=0, 2/2(0) = 1/20, Z2(0) = 0. ^ ' 

Subsequently, a xz- symmetric periodic orbit can be represented by a point 
in the four- dimensional space of initial conditions 

n4 = {(xio, a;2o, 2^20, 2/20)} • 

If an orbit has two perpendicular crossings with the then it is 

symmetric with respect to that axis. We can assume that at t = both 
planets start perpendicularly from the x-axis and the periodic orbit has initial 
conditions 

a;i(0) = xio, ^2(0) = X20, 2/2(0) = 0, ^2(0) = 0, , . 

i;i(0) = 0, X2(0)=0, 1/2(0) = 2/20, i2(0) = i20. ^ ' 

Thus, a x-symmetric periodic orbit can be represented by a point in the 
four-dimensional space of initial conditions 

n4 = {(xio, X20, 2/20, 2^20)}- 

By changing the value of Z2 (for the x2;-symmetry) or Z2 (for the x-symmetry) 
a monoparametric family of periodic orbits is formed. Also, a monoparamet- 
ric family can be formed by changing the mass of a planet. Pi or P2, but 
keeping the value Z2 (or Z2) constant. 



2.3 Linear stability of periodic orbits 

The linear stability of the periodic orbits can be found by computing the 
monodromy matrix M of the variational equations of the system Eqs. (|6|- 
([9]) . Since the system is of four degrees of freedom, M is an 8 x 8 symplectic 
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matrix and has four pairs of conjugate eigenvalues. One pair is always equal 
to unity, because of the existence of the energy integral. In order for the 
periodic orbit to be linearly stable, all of the eigenvalues have to be on the 
unit circle. Single instability is considered if there exists one pair of real 
eigenvalues {/i, If an eigenvalue /i is complex, then there exist also 
the eigenvalues /i, fi~^ and and the stability is characterized as complex 
instability, while the remaining eigenvalue pair may either lie or not on the 
unit circle ( [Marchal 1990; Hadjidemetriou 2006b). Nevertheless, apart from 



the pair (1,1), M may have two or three pairs of real eigenvalues. 
2.4 Periodic orbits and resonances 

When the planetary masses are very small compared to the mass of the Star, 
periodic orbits in the rotating frame correspond to almost Keplerian plane- 
tary orbits in space with orbital elements (semimajor axis), (eccentric- 
ity), ii (inclination), Ui (argument of pericenter), f2j (longitude of ascending 
node) and Aj (mean longitude) or Mj (mean anomaly), where i = 1,2 indi- 
cates the particular planet. 

A periodic orbit defines an exact resonance for the underlying dynamical 
system and is locked to a mean motion resonance np/nj = {p + q)/p-, where 
nj is the mean motion of the Jupiter and np is the mean motion of the other 
planet. In this paper, we study the 1/2 (Jupiter is the inner planet Pi) or 
the 2/1 (Jupiter is the outer planet P2) resonant cases. We distinguish 2/1 
from 1/2 resonance because different dynamics is obtained for the restricted 
problem. However, this distinction is not essential for the general problem. 
We can define the resonant angles ai and a2 of the system as: 

qai = {p + q)\i- p\2 - q^i 
qa2 = (p + g)Ai - p\2 - qTU2 

In the following, we refer to the resonant angles cti, Aw=(Ti — (T2 and AQ = 



^2 — ^^1- For the initial conditions (11) the orbital elements Ui and Qi can 
take the values | or For the initial conditions (12) the orbital elements 



uji and fli can take the values or vr. Also, for both symmetries the planets 
at t = are at their periastron or apoastron. Subsequently, the resonant 
angles that correspond to a periodic orbit can take values or tt. 
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3 Resonant Periodic Orbits and planetary evo- 
lution in their vicinity 

The computation of periodic orbits is based on the satisfaction of the peri- 
odicity conditions 

xi(T/2) = 0, 0:2(^/2) = 0, Z2{T/2) = {xz- symmetry) 

or 

Xi{T/2) = 0, X2(r/2) = 0, Z2iT/2) = (x - symmetry) 

where T/2 is the time of the first crossing of the section plane j/2 = 0. Keeping 
Xio fixed in the set of initial conditions, 114 or 114, we differentially correct 
the rest initial conditions, until they satisfy the periodicity conditions with 
accuracy of 11-12 digits. 

As already stated, the stability of a periodic orbit is defined by the 3 
pairs of non unit eigenvalues of the monodromy matrix. However, due to 
the limited accuracy, we cannot be sure, whether the eigenvalues lie on the 
unit circle, or not. An estimation of the accuracy of computations can be 
retrieved by the accuracy of computation of the unit eigenvalues which is 
about 10 digits. Nevertheless, apart from the linear stability, we estimate 
the evolution stability by using as index the Fast Lyapunov Indicator (FLI), 



(Froeschle et al. 1997) and, particularly, we compute the de-trended FLI, 



(Voyatzis 2008), defined as 



DFLI{t) = log Qmax{|ei(t)|, 

where are deviation vectors (initially orthogonal) computed after numerical 
integration of the variational equations. 



3.1 Examples of periodic orbits 

In Fig. [l| we illustrate an 1/2 resonant stable xz-symmetric periodic orbitj^in 
the inertial frame and in the projection space X2y2Z2 of the rotating frame and 
for planetary masses mi = 10~^ and m2 = 4 x 10~^. The initial conditions 
correspond to the planetary orbital elements 

ai = 0.45, ei = 0.66, ii = 12.79°, = 270°, coi = 90°, Mi = 0° 
a2 = 0.71, 62 = 0.34, i2 = 20.46°, ^2 = 90°, 102 = 270°, M2 = 180°. 



^This orbit belongs to the family Fg(1, see Section 



4.2 
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Its xz symmetry is revealed by the projections of the orbit in the planes of 
the space X2y2Z2 (see Fig. [T]d). The corresponding resonant angles take the 
values ai = 0°, Atu = 0° and Afi = 180°. 

In Fig. [2| we show an 1/2 resonant x-symmetric unstable periodic orbilj^ 
The planetary masses are mi = 10~^ and m2 = 10~^ and the initial conditions 
correspond to the planetary orbital elements 



ai = 0.33, ei = 0.42, ii = 0.26°, fii = 0°, ui = 0°, Mi = 0° 
02 = 0.52, 62 = 0.20, i2 = 19.79°, ^2 = 0°, uj2 = 0°, M2 = 180°. 



Its X symmetry is revealed by the projections of the orbit in the planes of 
the space X2y2Z2 (see Fig. ^p). The corresponding resonant angles take the 
values ai = 0°, Aro = 0° and A(] = 0°. 

The distribution of eigenvalues of the linear stability analysis is presented 
in Fig. [3] Regarding the stable a;2;-symmetric periodic orbit, all eigenvalues 
lie on the unit circle (Fig. |3^). One pair (of non unit eigenvalues) is very 
close to 1 and this is shown in the corresponding magnification (Fig. 
Concerning the unstable x-symmetric orbit there exists one pair of real eigen- 
values indicating the instability (Fig. [3|d). Also, one pair of eigenvalues is 
very close to 1 (see the magnification in Fig. |3ji). 

In Figjlj we present two different distributions of ei genv alues obtained for 



1/2 

two orbits that belong to the family Fg2i (see Section 4.2). In the first case, 
we have the presence of complex instability, while in the second one, two pairs 
of real eigenvalues exist. Along a family of periodic orbits the eigenvalues 
move smoothly either on the unit circle, or the real axis and can either enter 
the real axis, or the unit circle. Therefore, the stability type can change 
along the families. In the rest of the paper, we will indicate whether an orbit 
is stable or unstable without declaring the particular type of instability. 



3.2 Stability of evolution in the vicinity of periodic 
orbits 

It is known, that in the phase space, a stable periodic orbit is surrounded by 
invariant tori, where the motion is regular and stays in the neighbourhood of 
the periodic orbit. In the case of an unstable periodic orbit, chaotic domains 
exist at least close to the periodic orbit. In general, starting from the neig- 
bourhood of an unstable periodic orbit, the evolution of planetary system 
shows instabilities and chaotic behaviour. Numerical integrations show that 
such instabilities appear stronger as far as highly eccentric orbits are con- 



^This orbit belongs to the family G^(^j, see Section 



4.2 
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(a) (b) 

Figure 1: A symmetric periodic orbit with respect to the a;z-plane. The 
planetary mass ratio p = 1112/1111 equals to 0.4. a The planetary orbits in the 
inertial frame (the mutual inclination is 33°). b The periodic orbit in 

the rotating frame together with its projection to the planes of the frame. 




Figure 2: A symmetric periodic orbit with respect to the x-axis. The plan- 
etary mass ratio p = 1712/1711 equals to 0.01. a The periodic orbits of the 
planets in the inertial frame (the mutual inclination is Ai ^ 20°). b The 
periodic orbit of the third body together with its projection to the planes in 
the rotating frame. 
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Figure 3: Presentation of the linear stability, a The eigenvalues of the stable 
periodic orbit of Fig. [1} b The eigenvalues of the unstable periodic orbit of 
Fig. |2} c and d Magnifications of cases (a) and (b), respectively. 





Figure 4: An example of transition from a complex instability to b real 

1 /2 

instability along the family F^g.i P = 0-4 (see Section 4.2). 
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cerned (Voyatzis 2008). Hereafter, we present some numerical results that 
show the above mentioned behaviour of orbits. 

We consider the initial conditions of the unstable x-symmetric periodic 
orbit given in the previous section (Fig. [2]). Since the initial conditions of 
the periodic orbit are not exact, the instability leads the evolution far from 
the periodic orbit. In Fig. [5} we present the evolution of the orbital elements 
Oj, Cj and and the resonant angles ai, Aw and AQ. We observe that 
the eccentricity and the inclination of the heavier planet Pi remain almost 
constant during the whole time of integration. However, the planet P2 after 
the critical time 12 x 10^ t.u. destabilizes, its eccentricity increases and then 
oscillates around a high value. Also, its inclination takes large values and 
shows large oscillation around 90°, i.e. its motion turns from prograde to 
retrograde and vice versa. A significant oscillation of the semimajor axis of 
P2 after 12 x 10^ t.u. is also, apparent, though the system remains in the mean 
motion resonance. The resonant angles cxi and Aw initially librate arround 
0°, while AQ increases. After the critical time, ai and Aw start to rotate, 
while AQ librates irregularly in the domain [52°, 145°] . The computation 
of the DFLI (see Fig. [sj case II) shows clearly the chaotic nature of the 
evolution. 

If we start with the initial conditions given for the stable xz-symmetric 
periodic orbit (Fig. [T]), we will obtain that the orbital elements Oj, e^, ii and 
the resonant angles remain constant for ever (assuming sufficient integration 
accuracy). Starting near the periodic orbit, and particularly considering a 
deviation of 5° in Qi (i.e. we set now Qi = 275°), we obtain a regular 
evolution as it is shown in Fig. |6| The inclination and the eccentricity values 
of both planets oscillate arround their initial values, while the resonant angles 
cTi, Aw and AQ show librations around the values which correspond to the 
periodic orbit, namely 0°, 0° and 180°, respectively. Now, the DFLI does not 
increase (Fig. [sj case I) 

If we start not sufficiently close to the stable periodic orbit, the stability 
of evolution is not guaranteed. We consider a larger deviation, compared to 
the deviation used in the previous case. Particularly, we consider a deviation 
of 15° in ^1 and we integrate the orbit. The results are shown in Fig. [7} 
The chaotic behaviour of the evolution is obvious and finally the planet P2 is 
scattered. The resonant angles ai and Aw rotate, but it is remarkable that 
AQ librates irregularly around 180°. The DFLI clearly shows the chaotic 
evolution (Fig. Isl case HI). 
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20000 
t 



(b) 



Figure 5: Evolution of orbital elements and resonant angles starting with 
initial conditions (of limited accuracy) of the unstable periodic orbit of Fig. 
|2} Red and black line stands for Jupiter (Pi) and planet P2, respectively. 




40000 80000 120000 40000 80000 120000 

t t 

(a) (b) 

Figure 6: Evolution of orbital elements and resonant angles of an orbit which 
starts in the vicinity of the stable periodic orbit of Fig. [TJ Red and black 
line stands for Jupiter (Pi) and planet P2, respectively. 
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Figure 7: Evolution of orbital elements and resonant angles of an orbit that 
starts relatively far from the stable periodic orbit of Fig. [Tj The orbit is 
chaotic and after a time span close encounters between planets cause the 
scattering of the outer planet P2- Red and black line stands for Jupiter (Pi) 
and planet, P2, respectively. 



100 




400 800 1200 1600 2000 

t 



Figure 8: Computation of the DFLI. Case I corresponds to the regular evo- 
lution of Fig. |6} case II and case III to the chaotic evolution presented in 
Figs. |5]and[7| respectively. 
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Figure 9: Symmetric periodic orbits of the 3D-CRTBP in a 2:1 resonance 
and b 1:2 resonance. Blue color corresponds to stable periodic orbits and 
the red one to unstable orbits. The initial conditions of orbits presented by 
dots used as starting points for the continuation in the general problem. 



4 Families of periodic orbits 

We compute families of periodic orbits for the general spatial problem (3D- 
GTBP) by analytical continuation of periodic orbits either from the spatial 
circular restricted problem (3D-CRTBP), called Scheme I, or from the pla- 
nar general problem (2D-GTBP), called Scheme II. The presentation of such 
families will be given in projection planes defined by the planetary eccentric- 
ity and inclination values. 



4.1 Scheme I: Continuation from 3D-CRTBP to 3D- 
GTBP 

We consider the restricted problem with Jupiter moving in a circular planar 
motion with period 27r and a massless planet which moves in the space. 
In this problem, we can obtain families of periodic orbits which bifurcate 



from vertical critical periodic orbits of the planar restricted problem (Henon 



1973). These families exist either when Jupiter is the inner, or the outer 



planet (resonance either 1/2, or 2/1, respectively). 



Following the notation of Kotoulas and Voyatzis (2005), families of pe- 



riodic orbits which are symmetric with respect to the x^-plane are denoted 
by F and families of periodic orbits symmetric with respect to the 
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are denoted by G. The subscript c or e indicates the spatial restricted prob- 
lem (circular or elliptic, respectively) from which these families bifurcate. 
Moreover, m indicates the continuation with respect to the mass of the ini- 
tially massless body. Also, the corresponding resonance is indicated by a 
superscript. 

We computed and present the families of 2/1 and 1/2 resonant symmet- 



ric periodic orbits of the 3D-RTBP in Fig. 9 (see also Hadjidemetriou and 



Voyatzis (2000) and Kotoulas (2005), respectively). We denote these fami- 
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lies by F ^ , where the subscript i indicates the continuation from plane to 



space (see Section 4.2). For the 2/1 resonance such families correspond to 
relatively high eccentricity values of the massless planet Pi. Conversely, in 
the 1/2 resonance all periodic orbits correspond to almost circular orbits of 
the massless planet P2 and are unstable. 

It is known that the periodic orbits of the 3D-RCTBP of period T can 
be analytically continued, by increasing the planetary mass of the initially 
massless body to a periodic orbit of the 3D-GTBP of the same symmetry. 



provided that T 7^ 2k-K, where k is an integer (Ichtiaroglou and Michalodimi- 
trakii|[l980| . This is the case for all orbits of the families given in Fig. |9| For 



the particular computations presented hereafter, we continue analytically the 
periodic orbits indicated by the dots in Fig. [9] with respect to the mass mi 
(inner planet) or m2 (outer planet) according to the resonance 2/1 or 1/2, re- 
spectively. These periodic orbits correspond to some initial inclination value 
of the massless body, iio or 220, and identify the generated families. 

2/1 2/1 

4.1.1 Families Fc'm- Continuation from orbits of family F^ - with 
respect to mi 



In Fig. 10, we present the continuation of the x^-symmetric periodic orbits 



of family F^{^. The initial inclination values, iio, are indicated. The formed 

2/1 

families are denoted by Fc'm and their periodic orbits correspond to the 
resonant angles 

Aw = 0°, An = 180°, ai = 0°. 

The variation of the eccentricity of the planet Pi along the families with 
respect to the mass mi is shown in Fig. 10 1. As mi increases from to 10~^ 



(approximately), the eccentricity, 62, of Jupiter increases rapidly from to 
0.25 for all computed families. Then, for small initial inclinations, iiQ, of Pi, 
62 remains almost constant along the families, but continues to increase for 
higher initial inclinations, zio (see panel b). The diagram rrii-ii (panel c) 
shows a plateau for the inclination of the inner planet, whose range extends 
up to almost mi = 0.01 for the families of small inclination values of Pi. In 
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Figure 10: Families Fc'm of periodic orbits which are continued from the 
restricted problem (3D-CRTBP) and are symmetric with respect to xz-plane. 
The families are presented by characteristic curves in different projection 
planes. For each family the inclination value, zio, of the starting orbit is 
indicated. Bold blue and red coloured segments represent stable or unstable 
orbits, respectively. 
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Figure 11: The variation of mutual inclination of planets along the stable 
segments of families Fclm- The parameter value iio is shown for each family. 

contrary, Jupiter (P2) shows an almost linear increment of its inclination, as 

mi increases (panel d). 

2/1 

Families Fc'm start having stable periodic orbits. The stability extends up 
to critical value of mi, which depends on the particular initial value iiQ. In 
these stable segments the mutual planetary inclination varies, as it is shown 



in Fig. [TT| We observe that we can have stable planetary orbits with mutual 
inclination, Ai, up to approximately 50° - 60° when the inner planet becomes 
very massive with respect to Jupiter. 

4.1.2 Families Fljm • Continuation from orbits of family fI^^ with 
respect to m2 

The continuation of the 1/2 resonant x2;-symmetric periodic orbits that be- 

1 /2 

long to the family Fj^ of the restricted problem is presented (in a similar 

way as above) in Fig. x.^ uiul-c tuc laiimj 10 unouauic, tuc laiiimco -t c,m, 

lie analytical continuation with respect to m2, start 



Since the family F^^^ is unstable, the families Flf'^ 

which are obtained by t 
having unstable periodic orbits, too. However, as m2 increases the orbits 
with < 35° become stable. Instability becomes again apparent after a 
larger value of m2, which depends on the initial inclination, i2o- It is inter- 
esting to note that the characteristic curves, which represent the families in 
the projection planes, show a folding. Namely, for a particular value of m2 
we may obtain more than one periodic orbits that belong in the same family. 



This folding becomes clear in the magnified plot in each panel of Fig. 12 



The mutual planetary inclination along the stable segments is shown in Fig. 
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1/2 

Figure 12: Families Fc'm of periodic orbits presented in a similar manner 



as in Fig. 10 Each family is identified by the inclination value, ^205 of the 



starting orbit. 
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Figure 13: The variation of mutual inclination of planets along the stable 
segments of families Fljm- The parameter value ^20 is shown for each family. 



13| As in the previous case, we obtain stable planetary orbits up to relatively 

high mutual inclination values (approximately 50°). 

1/2 

As we have mentioned, the families Fc'm start with orbits where the ec- 
centricity of the planet P2 is very small. For 120 < 26° the corresponding 

1 /2 

resonant angles of the family [ of the restricted problem are 

(i) Aro = 180°, /\VL = 180°, ai = 0°. 

For > 26° the orbit of P2 changes apsidal orientation and the resonant 
angles take the values 

(ii) Aw = 0°, An = 180°, ai = 180°. 

1/2 

The orbits of the generated families Fc/m start with the above resonant angle 
values, (i) or (ii), accordingly. However, the families that start with the values 
(i), show a decrement of the eccentricity 62, as m2 increases. At 62 = we 
obtain a change of the apsidal orientation of P2 and the rest of the family 
consists of orbits with the resonant angle values (ii). 

2/1 2/1 

4.1.3 Families Gc'm- Continuation from orbits of family G^ j with 
respect to mi 



In Fig. 

erated 



14 



2/1 

we present the families Gc'm of x-symmetric periodic orbits gen- 
rom orbits of the unstable family of the restricted problem. The 
continuation of these periodic orbits, which is performed by increasing the 
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mass mi, extends only up to a certain value m*. Note that the value of ml for 
all families is relatively very small (of order 10~^) with respect to Jupiter's 
mass 777,2 = 10"'^. Starting from mi = m^ and by decreasing mi we can ob- 
tain a different branch of the family down to mi = 0. However, at this point, 
as far as Jupiter is concerned, we have 62 7^ and 72 = 0, i.e. the families 
end at periodic orbits of the elliptic restricted problem (3D-ERTBP). Note, 
now, that Jupiter's orbit is almost planar in all cases. 

2/1 

All orbits of families Gc'm are unstable. Their resonant angles are 
Azu = 0°, An = 0°, ai = 180°. 



1 /2 1 /2 

4.1.4 Families Gc'm '• Continuation from orbits of family G^ - with 
respect to m2 

1 /2 

Similar behaviour as above (start and end) is obtained for the families Gc'm, 



where now Jupiter is the inner planet (-Pi). AH orbits are unstable (Fig. 15). 

1 /2 

Concerning the resonant angles, the generating family • of the re- 
stricted problem starts with orbits having 

(i) Aw = 180°, An = 0°, (Ti = 180°, 

until 720 < 50°. For larger inclination values we have a change in the apsidal 
orientation and we get the values 

(ii) Aw = 0°, An = 0°, (Xi = 0° 

1 /2 

The generated families Gc'm start with the values as above, (i) or (ii) for 
720 < 50° or 720 > 50°, respectively. For case (i), the eccentricity 62 starts 
decreasing down to zero along the family. Then, the apsidal orientation of 
planet P2 changes and the resonant angles take the values (ii). 



4.2 Scheme II: 
GTBP 



Continuation from 2D-GTBP to 3D- 



Resonant symmetric periodic orbits for the planar general problem (2D- 
GTBP) can be obtained by starting from either the circular family of orbits 



(e.g. Hadjidemetriou (2006a)), or from resonant families of the planar re- 
stricted problem (e.g. Antoniadou et al. (2011)). Particularly, for the 2/1 



(or 1/2) resonance, planar periodic orbits have been computed by 1 


Beauge 


et al. 2006 


Voyatzis and Hadjidemetriou 2005 


Voyatzis et al. 


2009 


). Now, 



the distinction between the 2/1 and 1/2 resonance is not essential, since both 
planetary masses are nonzero and we refer to the mass ratio p = 1712 /mi, 
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where Pi is the inner planet. In the following computations, we always set 



mi = 0.001 (Jupiter's mass). In Fig. 16, we present the 1/2 resonant families 



fi and /2 of the general planar problem. The families are given for various 
values of the mass ratio p in the projection space defined by the planetary 
eccentricities. 

The variational equation that corresponds to the z-component of the 
motion of planet P2 is written a^ 

(2 = ^(2 + 5(2 (13) 

where 

mmo[l- "'2('^^2+!/2) ] mmi[l+ '"°'"2(9^2+^2) ] 

m6x-[ mm-\ 6x-\ 



(14) 

If A(T) = {^ij}, i,j = 1,2, is the monodromy matrix of Eq. 13 for a 



TJ _ rnom2y2 mom2y2 

^ (rrio+mi)toi[(xi-X2)2+J/2]3/2 [(!r^l^+^2)2+y2]3/2 • 



planar periodic orbit of period T, then, this orbit has a vertical critical index 

(15) 



( Michalodimitrakis 1979a) equal to 

^ll'C22 + ^21^12 



^1^22 — ^21^2 

Any periodic orbit of the planar problem with \ay\ = 1 is vertical critical 
(v. CO.) and can be used as a starting orbit for analytical continuation in the 
spatial problem with the inclination of P2 as a parameter (in computations 
we increase Z2 or ^2)- 

In Fig. 16, we present the v.c.o. by dots (see also Tables fT]|3]) . In Fig. 



rffi , we illustrate the variation of the index ay along the family /i for some 
values of p. For p < 0.0372, we find three v.c.o. and starting from one of 

them (see Table ITl) we can analytically continue x^-symmetric periodic orbits 

1/2 

in space. We denote such families by F^/ • where the subscript g indicates 
that the family bifurcates from the general planar problem (it is followed by 
a number that states the generating planar family) and i indicates that the 
continuation takes place from plane to space. From the remaining couple 
of v.c.o. (see Table [2]) x-symmetric periodic orbits in space are obtained 
forming the families G]^{\. An example of the above families is given in Fig. 



18 1 in the projection space 6162 Ai, where Ai is the mutual inclination. We 



1 /2 

obtain that the family Gg[ ^ starts from and ends at v.c.o. of the same planar 
family /i. For p > 0.0372 the couple of v.c.o., from which we obtained the 
family G]^{\ disappears. This transition is explained by the variation of the 
index ay along the family /i and is shown in Fig. [17) 3. 
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(a) 



(b) 



Figure 16: Families of symmetric periodic orbits of the planar general prob- 
lem a fi and b /2. Bold blue colour stands for stable periodic orbits and red 
for the unstable ones. 




Figure 17: a The variation of of the index a„ along the family /i (ei is used 
as the parameter of the family). For p < 0.0372 the index becomes critical in 
three points, while for p > 0.0372 we obtain only one vertical critical point, 
b A magnification of panel (a) which illustrates the disappearance of the 
couple of critical orbits. 
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p=0.02 



p=OM 




p=0.1 



p=OA 



Figure 18: Bifurcation of families F^lJ (and G^i^j only for p < 0.0372) which 



are generated from vertical critical points of families /i given in Fig. 16 1. 
Bold blue colour stands for stable periodic orbits and red for the unstable 
ones. 
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Table 1: Eccentricity values for the v.c.o. of families /i from which xz- 
symmetric periodic orbits bifurcate in space. 



p 


ei 


62 


0.01 


0.014402 


0.047445 


0.1 


0.442008 


0.198520 


0.2 


0.5486138 


0.2496163 


0.3 


0.5750049 


0.2596675 


0.4 


0.590191 


0.264541 


0.6 


0.6051763 


0.2672597 


0.8 


0.6131065 


0.2670630 


1 


0.6182749 


0.2659223 


2 


0.6292770 


0.2563407 


5 


0.6350957 


0.2287521 


10 


0.6376757 


0.1995789 


20 


0.6375954 


0.1613843 



Table 2: Eccentricity values for the v.c.o. of families /i from which x- 
symmetric periodic orbits bifurcate in space. 



p 


ei 


62 


0.01 


0.046879 
0.451260 


0.030487 
0.208753 


0.02 


0.099993 
0.420425 


0.004378 
0.191886 


0.03 


0.176672 
0.366875 


0.052130 
0.162298 



1/2 

Concerning the families Fg{^, up to the mass ratio value p* ~ 0.12 all 
the corresponding v.c.o. are unstable and, subsequently, the continuation 
generates families that start as unstable. The families continue to be unstable 
until the end of the computations. For p > p* the v.c.o. belong to the stable 

1 /2 

segment of family /i. Now the generated families Fg{- start as stable and 
show a monotonical increase of the inclinations ii and i2, and, subsequently. 



Ai, too (see Fig. 18). The resonant angles that correspond to these periodic 
orbits are 

Avj = 0°, An = 180°, di = 0°. 

1 /2 

The stability type along the families Fg{ ■ changes above a maximum value 



•^Notc that the variational equation given in ( Michalodimitrakisl|l979a ) holds only for 
the normalization mo + mi = 1. 
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Table 3: Eccentricity values for the v.c.o. of families /2 from which xz- 
symmetric periodic orbits bifurcate in space. 



p 


ei 


62 


0.01 


0.003953 


0.461514 


0.1 


0.037795 


0.447835 


0.2 


0.072314 


0.435873 


0.3 


0.104294 


0.426776 


0.4 


0.133030 


0.414625 


0.6 


0.185104 


0.397536 


0.8 


0.231733 


0.386791 


1 


0.272251 


0.375495 


2 


0.431450 


0.354297 


5 


0.683522 


0.374495 


10 


0.813132 


0.419156 


20 


0.880779 


0.437977 



for the mutual inclination, Ai^ax that depends on p, as it is shown in Fig. 
[19} In case (a), Jupiter is the inner planet, Pi, and P2 takes large mass values 
(up to 20 Jupiter masses), as p increases. In case (b), we present the same 
computations, but now Jupiter is the outer planet, P2 and the other planet. 
Pi, takes smaller and smaller mass' values, as p increases. We observe that 
we can have stable orbits up to mutual inclinations 40° < Aimax < 50° and 
there is no essential difference, whether Jupiter is the inner, or the outer 
planet in the system. 



1/2 

In Fig. 20, we present the families Fg2i, which bifurcate from the v.c.o. 



of families jTgiven in Fig. 16 d. The corresponding resonant angles are 



Azu = 180°, An = 180°, ai = 180°. 

The families which start from stable v.c.o. (in the range 0.3 < p < 7, 
approximately) show a segment of stable periodic orbits up to a critical 
mutual inclination value, Aimax- Aimax seems to increase, as p increases in 
the above mentioned interval and reaches 30°, approximately. The unstable 
segments show transitions from real to complex instability and vice versa 
(see e.g. Fig. |4]). 

5 Miscellaneous cases 



2/1 

The families Fc'm are presented in Fig. 10 only for mi < 0.01. Further 



continuation with respect to the mass, mi, is possible and as a result, we 
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|30- 



(a) 



(b) 



(a) mi=10"^, in2=variable 

(b) m\=variable, m-f=lO~^ 



Figure 19: The maximum mutual inclination Ai^ax as a function of mass 
ratio p = 1712 /nil up to which we obtain stable periodic orbits in families 

1 /2 

j. Stability exists for p > 0.12. In case (a) is mi = 0.001 (Jupiter is the 
inner planet, Pi) and in case (b) m2 = 0.001 (Jupiter is the outer planet, P2) 




Figure 20: The families Fg2i which bifurcate from the v.c.o. of families /2 of 
the planar general problem. 
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Figure 21: a Families Fcjm up to large mass values of Pi b Families Gc,m 

2/1 

and Ge\m computed by analytical continuation with respect to mi. 



obtain the characteristic curves of Fig. 21 1. We observe that for mi > 0.01 

* * 1 /2 

foldings become apparent as in families Fc'm (see Fig. 12|). Also, we can 



obtain the formation of a loop, given by the separatrix characteristic curve, 
which includes closed families. Thus, along a family Fc/m we may obtain 
more than one periodic orbits (e.g. the periodic orbits A, B and C) that 
correspond to the same planetary masses. 

2/1 

We have computed the families Gc'm starting from generating orbits of 
the 3D-CRTBP. This explains the subscript c in their name, although these 
families end at orbits of the 3D-ERTBP. All of these orbits extend up to a 
critical mass value, ml and correspond to inclinations iio ^ 70°. However, if 
we consider as starting points periodic orbits of the 3D-ERTBP with iiQ > 70° 
then, by continuation with respect to the mass, mi, we obtain the families 

2 /I 

Ge,m which do not end at the circular restricted problem. This is illustrated 



by the characteristic curves presented in Fig. pTp. We found that a similar 
situation holds, also, for the 1/2 resonance, namely beside the families Gc'm 
(see Fig. 15| there exist, also, families Gefm- All of them are unstable. 

The families obtained by the analytical continuation of Scheme I, consist 
of periodic orbits with Z2 ^ or Z2 ^ 0, namely the planetary orbits are of 
nonzero inclination. Assuming any of the above mentioned periodic orbits we 
can apply the continuation of Scheme II, in order to decrease the inclination 
and terminate in the planar general problem. Then, this terminating point 



must be a v.c.o. of the planar problem. An example is shown in Fig. 22 
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1 /2 

Figure 22: An example of a family Fg{^ which appears as a bridge between 
a planar family fi and a family Fc/m. Also, a case of a family F^[^, which 

2/1 

starts and ends at orbits of family Fc/m, is shown. 

2/1 

We start from the periodic orbit C that belongs to the family Fc'm for mi = 

2/1 

0.0139. Along Fc,m it is Z2 = 0.01042 (constant). By decreasing Z2, according 
to the Scheme II, we reach the planar general problem and a v.c.o that belongs 

1 /2 

to the family /i with p = 0.07. Thus, we obtain a F^/j family. However, 
this is not always the case. Considering, for example, as a starting periodic 
orbit, the orbit A, the continuation by varying Z2 and keeping mi constant 
results to a family that terminates at the orbit B without reaching the plane 
Z2 = 0. This family evolves exclusively in the 3D-GTBP and we denote it 

2/1 

by F- [ . The existence of a folding (namely two points corresponding to the 
same planetary masses) is a necessary condition for obtaining such a type of 
families. 

6 Conclusions 

In this paper, we presented 2:1 resonant families of periodic orbits in the three 
dimensional general TBP of planetary type considering a suitable rotating 
frame of reference, which reduces the system to four degrees of freedom. 

Symmetric periodic orbits are of two types: symmetric with respect to 
the xz-'plane of the rotating frame (Afi = 180°) and symmetric with respect 
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to the X-axis {AQ = 0°). For the computation of famihes of periodic orbits 
we can apply two continuation schemes: 

i) we start from periodic orbits of the 3D circular or elliptic restricted TBP 
and continue them by varying the mass of the small planet (which is initially 
the massless body) 

ii) we start from vertical critical periodic orbits that belong to the families 
of the general planar problem and continue them in space by varying the 
planetary inclination. 

The computation and continuation procedures can be applied similarly 
to any other resonance. Moreover, similar procedures can be used for the 
computation of asymmetric periodic orbits too, but in this case, it is required 
to meet a larger number of periodic conditions. Also, we considered and 
presented orbits only for prograde planetary motion. 

Most of the periodic orbits found are linearly unstable. Especially, all 
periodic orbits symmetric with respect to the of x-axis are unstable. In their 
neigbourhood in phase space, we have planetary orbits that evolve chaoti- 
cally. We obtain relatively large excitations in eccentricity or inclination and 
after a relatively short-term evolution the planetary system is disrupted, due 
to close encounters. 

Nevertheless, families of a;2;-plane symmetry show segments which include 
stable periodic orbits. Such orbits can be obtained up to 50° — 60° of mutual 
planetary inclination. Particularly we found that all families which bifurcate 
from stable vertical orbits of the planar problem start having stable periodic 
orbits. For the case of the planar family fi (cri = 0°, Aw = 0°), stability 
of the bifurcated 3D periodic orbits is observed for a wide range of values of 
the planetary mass ratio, particularly for p > 0.12 and at least up to p = 20 
(where we stopped our computations) . For the families of 3D periodic orbits 
which bifurcate from the planar family /2 (o"i = 180°, Accj = 180°) stable 
orbits exist for 0.3 < p < 7. The maximum mutual inclination which is 
observed for stable orbits increases as p increases in the above domain and 
reaches 30°, approximately. All stable periodic orbits found correspond to 
quite eccentric planetary orbits (for at least one of the planets) . 

Starting with initial conditions sufficiently close to a stable periodic orbit 
the planetary evolution takes place regularly. The planetary system remains 
in the mean motion resonance and the eccentricities and inclinations show 
small oscillations. Also, all resonant angles hbrate around the values that 
correspond to the stable periodic orbit. Such regions near stable periodic 
orbits may be considered as strong candidates for hosting planetary systems. 
The relation between the dynamics of real planetary systems and 3D periodic 
orbits of the general TBP remains to be studied. 
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